################################################################
# Replication File: "The Politics of Too Much - The Decline in Labor Power 
# and the Global Rise of Corporate Saving" 
# Author: Nils Redeker
# Firm level analysis replication
# R version 3.6.2 (2019-12-12)
#######################################################################

# clean working environment
rm(list = ls())

# load packages --- --- ---- --- --- --- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- --- --- --- --- 

library(readr)
library(rdd)
library(rdrobust)
library(ggplot2)
library(stargazer)
library(tidyr)
library(dplyr)
library(RPostgres)
library(ggpmisc)
library(broom)
library(lfe)
library(readxl)

# load data --- --- ---- --- --- --- --- --- --- --- --- --- -- --- --- --- --- --- --- -- --- --- --- --- ----- 

# set working directory
setwd()

# Firm-level Compustat data is available from WRDS. For instructions on how to 
# acquire the data, please refer to README file. 

firms.de<- read_csv("~/JOP_replication_firm_level.R") 

#------------------------------------------------------------------------------------------------------#
# I) MAIN ANALYSIS
#------------------------------------------------------------------------------------------------------#

#------------------------------------------------------------------------------------------------------#
# 1) Prepare data for analysis
#------------------------------------------------------------------------------------------------------#

firms.de.analysis<-firms.de %>% filter(emp>1.5 & emp<2.5, fyear>1990 & fyear<2018) %>% # filter firms around threshold
  mutate(per=ifelse(emp>=2,1,0), # generate treatment dummy
         cash=ch + ivst, # generate cash variable 
         ln_cash=log(cash+1),
         cash_ratio = cash/at, # cash per total assets
         dvt_at = dvt/at, # dividiends 
         xstfws_at=xstfws/at, # staff expenses (wages & salaries)
         xstfo_at=xstfo/at, # staff expenses (other)
         lppent=dplyr::lag(ppent, n=1L),
         lintan=dplyr::lag(intan, n=1L),
         lat=dplyr::lag(at, n=1L),
         ldp=dplyr::lag(dp, n=1L),
         lcapx=dplyr::lag(capx, n=1L),
         capx.at=capx/lat,
         inv.grw=((capx-lcapx)+(intan-lintan))/lat, # investment changes (Peters & Taylor 2017)
         emp=emp*1000, 
         cash_emp=cash/emp, # cash per employee
         fyear=as.factor(fyear),
         sic_one=substr(sic,1,1),
         man=ifelse(sic_one==1|sic_one==2,1,0),
         service=ifelse(sic_one==7 | sic_one==8,1,0),
         tech=ifelse(sic_one==4 ,1,0),
         trade=ifelse(sic_one==5 ,1,0)) %>% # create one-digit SIC code
  filter(sic_one!=6, !grepl(' SE', conm), # drop financial firms 
         !gvkey %in% c("236260","242158","216075")) # drop SEs


#------------------------------------------------------------------------------------------------------#
# 2) FIGURE 3: REGRESSION DISCONTINUITY PLOTS (95% CONFIDENCE INTERVALS)
#------------------------------------------------------------------------------------------------------#

summary(rdbwselect(firms.de.analysis$cash, # calculate optimal bandwidth (BW) for RDD
                   firms.de.analysis$emp, c=2000)) 

out<-rdplot(firms.de.analysis$cash, firms.de.analysis$emp, # generate data for RDD plot
            c=2000, binselect = "esmv"
            ,x.lim=c(1800,2200), p=1, scale=2, shade=TRUE)

visuals_main<-data.frame(out$vars_bins) %>% 
  mutate(per=as.factor(ifelse(rdplot_mean_x>2000,1,0))) %>% # treatment dummy for plot
  filter(rdplot_mean_x>=1840 & rdplot_mean_x<=2160) # drop obs. outside optimal BW

## FIGURE 3a) linear fitting

ggplot(visuals_main, aes(x = rdplot_mean_x, y = rdplot_mean_y, colour=per)) + geom_point(alpha = 0.8, size=2.5)+
  scale_color_manual(values=c("#FE4365", "#83AF9B")) + stat_smooth(method = "lm", level=.90, fill="LightCyan3", data=subset(visuals_main, visuals_main$per==1)) +
  stat_smooth(method = "lm", level=.95, fill="red", data=subset(visuals_main, visuals_main$per==0)) + ylim(-40, 150) +
  ylab("Corporate Savings (Million US$)") + xlab("Number of employees") +
  geom_vline(xintercept = 2000, linetype="dashed", 
             color = "black", size=0.5) + 
  annotate("text", x=1900, y=140, label= "No Parity Co-Determination", family="Frutiger-Light", size=6) +
  annotate("text", x=2100, y=140, label= "Parity Co-Determination", family="Frutiger-Light", size=6) +
  theme_linedraw() + theme(text=element_text(size=20, family="Frutiger-Light"), legend.position = "",
                           axis.text=element_text(size=20),
                           axis.title=element_text(size=20))

## FIGURE 3b) quadratic fitting

ggplot(visuals_main, aes(x = rdplot_mean_x, y = rdplot_mean_y, colour=per)) + geom_point(alpha = 0.8, size=2.5)+
  scale_color_manual(values=c("#FE4365", "#83AF9B")) + stat_smooth(method = "lm", formula = y ~ poly(x, 2), level=.95, fill="#83AF9B",
                               data=subset(visuals_main, visuals_main$per==1)) +
  stat_smooth(method = "lm", formula = y ~ poly(x, 2),level=.90, fill="#FE4365", 
              data=subset(visuals_main, visuals_main$per==0)) + ylim(-40, 150) +
  ylab("") + xlab("Number of employees") +
  geom_vline(xintercept = 2000, linetype="dashed", 
             color = "black", size=0.5) + 
  ylab("Corporate Savings (Million US$)") +
  annotate("text", x=1900, y=140, label= "No Parity Co-Determination", family="Lato-Light", size=6) +
  annotate("text", x=2100, y=140, label= "Parity Co-Determination", family="Lato-Light", size=6) +
  theme_linedraw() + theme(text=element_text(size=16, family="Lato-Light"), legend.position = "" ,
                           axis.text=element_text(size=16),
                           axis.title=element_text(size=16),
                           axis.text.y = element_blank())



#------------------------------------------------------------------------------------------------------#
# 3) TABLE 2: THE EFFECT OF LABOUR PARITY CO-DETERMINATION ON FIRM-LEVEL CASH-HOLDINGS
#------------------------------------------------------------------------------------------------------#
attach(firms.de.analysis)

# TABLE 2 - MODEL (1) - Baseline
model1<-(rdrobust(cash, emp, c=2000, all=TRUE))
model1$coef["Bias-Corrected",]
model1$pv["Bias-Corrected",]
model1$ci["Bias-Corrected",]


# TABLE 2 - MODEL (2) - Adding year/ indsutry-fixed effects
model2<-(rdrobust(cash, emp, c=2000, covs = cbind(fyear, man, service, tech, trade)))
model2$coef["Bias-Corrected",]
model2$pv["Bias-Corrected",]
model2$ci["Bias-Corrected",]

# TABLE 3 - MODEL (2) - Adding year/ indsutry-fixed effects
model3<-rdrobust(cash, emp, c=2000,
                  covs = cbind(fyear, man, service, tech, trade),
                  cluster=gvkey)
model3$coef["Bias-Corrected",]
model3$pv["Bias-Corrected",]
model3$ci["Bias-Corrected",]

#------------------------------------------------------------------------------------------------------#
# 4) FIGURE 4: PLACEBO TESTS WITH 95% CONFIDENCE INTERVALS
#------------------------------------------------------------------------------------------------------#


# Create data frame with placebo cut offs
placebo_cuts<-data.frame(rep(seq(1650,2350, by=25), each=3)) # generate alternative cut-offs
colnames(placebo_cuts)<-c("cuts")
true_cut<-data.frame(rep(2000, 3)) 
colnames(true_cut)<-c("cuts")
placebo_cuts<-rbind(placebo_cuts, true_cut)

# Prepare data for placebo RDDs
b <- seq(1650,2350, by=25) # generate alternative cut-offs
b<-append(b, 2000) # add true cut-off
container <- data.frame()
container2<-data.frame()

# Conduct placebo RDDs
for(i in b){
  a <- rdrobust(cash, emp, c=i)
  container <- rbind(container, a$coef)
  container2<-rbind(container2, a$se)
  placebo<-cbind(container, container2)
}

# Prepare data for plotting
placebo<-cbind(placebo, placebo_cuts)
colnames(placebo)<-c("coef","se","cuts")
ind <- seq(3, nrow(placebo), by=3)
placebo<-placebo[ind,] %>% mutate(upper95=coef+1.96*se,lower95=coef-1.96*se,
                                  real_cut=ifelse(cuts==2000,1,0))

# Display placebo tets
ggplot(placebo, aes(x = cuts, y = coef)) + 
  geom_point(data=subset(placebo, placebo$real_cut!=1),  shape = 21, fill = "grey", color = "grey",  size = 2) + 
  geom_errorbar(data=subset(placebo, placebo$real_cut!=1), aes(ymin=lower95, ymax=upper95), alpha=.4, fill="grey", width=0) +
  geom_hline(yintercept = 0, color= "black", linetype = "dashed") +
  geom_errorbar(data=subset(placebo, placebo$real_cut==1),aes(ymin = lower95, ymax = upper95), width = 15, color = "black") +
  geom_point(data=subset(placebo, placebo$real_cut==1), shape=21, color = "black", fill="black") +
  theme(text=element_text(size=12, family="Frutiger-Light"), legend.position = "") +
  xlab("Placebo Cutpoints") + ylab("Estimated Effect")


#------------------------------------------------------------------------------------------------------#
# 4) TABLE 3: MECHANISMS: EFFECT ON FIRMS’ SPENDING BEHAVIOUR
#------------------------------------------------------------------------------------------------------#

# Effect on dividends
model1<-rdrobust(dvt_at, emp, c=2000, 
                 cluster = gvkey,
                 covs = cbind(fyear, man, service, tech, trade))
model1$coef["Bias-Corrected",]
model1$pv["Bias-Corrected",]
model1$ci["Bias-Corrected",]
model1$bws

# Effect on staff expenses (wages & salaries)
model2<-rdrobust(xstfws_at, emp, c=2000, 
                 cluster = gvkey,
                 covs = cbind(fyear, man, service, tech, trade))
model2$coef["Bias-Corrected",]
model2$pv["Bias-Corrected",]
model2$ci["Bias-Corrected",]
model2$bws
# Effect on staff expenses (other)
model3<-rdrobust(xstfo_at, emp, c=2000, 
                 cluster = gvkey,
                 covs = cbind(fyear, man, service, tech, trade))
model3$coef["Bias-Corrected",]
model3$pv["Bias-Corrected",]
model3$ci["Bias-Corrected",]
model3$bws

# Effect on capital growth rate
model4<-rdrobust(inv.grw, emp, c=2000, 
                 cluster = gvkey, 
                 covs = cbind(fyear, man, tech, service, trade))
model4$coef["Bias-Corrected",]
model4$pv["Bias-Corrected",]
model4$ci["Bias-Corrected",]
model4$bws

#------------------------------------------------------------------------------------------------------#
# II) APPENDIX
#------------------------------------------------------------------------------------------------------#

#------------------------------------------------------------------------------------------------------#
# 1) FIGURE B5 - CASH HOLDING AND NET SAVINGS
#------------------------------------------------------------------------------------------------------#

# Create cash/ sic_one for all firms
firms.net_lend<-firms.de %>% mutate_all( ~replace(., is.na(.), 0)) %>% 
  mutate(fyear=as.numeric(as.character(fyear)), cash=ch + ivst, sic_one=substr(sic,1,1)) 

# Generating net lending variable based (see Appendix B1)
firms.net_lend<- firms.net_lend %>% mutate(opex=cogs + xsga,  
                                                 gos=sale-xopr+dp+xrd, 
                                                 gos=ifelse(is.na(gos)==0, ebitda + xrd, gos),
                                                 gs=gos - intpn - txt - dvt,
                                                 gfxf = aqc - sppiv + xrd,
                                                 net_lend = gs - gfxf)
# (a) 1990 - 2000: Cash holding and Net Savings
ggplot(firms.net_lend %>% filter(fyear>1990 & fyear < 2000, sic_one!="6"), aes(x= cash, y= net_lend)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", formula = y~x , se = F) +
  stat_poly_eq(aes(label = paste(..rr.label..)), 
               label.x.npc = 0.945, label.y.npc = 0.15,
               formula = formula, parse = TRUE, size = 7) +
  stat_fit_glance(method = 'lm',
                  method.args = list(formula = y~x ),
                  geom = 'text',
                  aes(label = paste("P-value = ", signif(..p.value.., digits = 4), sep = "")),
                  size = 7) + ylab("Net Lending") +
  xlab("Liquid Assets") + theme(text=element_text(size=16, family="Frutiger-Light"))

# (b) 2008 - 2015: Cash holding and Net Savings
ggplot(firms.net_lend %>% filter(fyear>2007 & fyear < 2016, sic_one!="6"), aes(x= cash, y= net_lend)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", formula =  y~x , se = F) +
  stat_poly_eq(aes(label = paste(..rr.label..)), 
               label.x.npc = 0.945, label.y.npc = 0.15,
               formula = formula, parse = TRUE, size = 7) +
  stat_fit_glance(method = 'lm',
                  method.args = list(formula =  y~x ),
                  geom = 'text',
                  aes(label = paste("P-value = ", signif(..p.value.., digits = 4), sep = "")),
                  label.x.npc = 0.95 , label.y.npc = 0.3, size = 7) + ylab("Net Lending") +
  xlab("Liquid Assets") + theme(text=element_text(size=16, family="Frutiger-Light"))  + ylab("") 

#------------------------------------------------------------------------------------------------------#
# 2) FIGURE B.6 - MCCRAY (2008) DENSITY TEST
#------------------------------------------------------------------------------------------------------#

DCdensity(emp, 2000, bin = NULL, plot=F)

ggplot(firms.de.analysis, aes(x=emp))+
  geom_bar(aes(x=emp,y=..density..), stat="bin", fill="grey") +
  theme(text=element_text(size=12, family="Frutiger-Light"), 
        axis.title=element_text(size=12)) + ylab("Density") +
  geom_vline(xintercept = 2000, linetype="dashed", 
             color = "black", size=0.7) +
  annotate("text", x=2300, y=0.0015, label= "P-value density test: 0.34", family="Frutiger-Light", size=4) +
  ylab("Density") +  xlab(expression(atop("Number of Employees", paste("   "))))


#------------------------------------------------------------------------------------------------------#
# 3) TABLE B.1 - FALSIFICATION TESTS: EFFECT OF PARITY CO-DETERMINATION ON PRE-TREATMENT COVARIATES
#------------------------------------------------------------------------------------------------------#

# load ownership data from ORBIS 
subset_ownership<- read_excel("~/Dropbox/02_PhD/07_Corporate_Savings/01Data/13Ownership_Orbis/Orbis_Ownership.xlsx")
subset_ownership<-subset_ownership%>% mutate_at(c("owner1_share","owner2_share",
                                                  "owner3_share", "owner4_share", 
                                                  "owner5_share"), as.numeric)

shares<-c("owner1_share","owner2_share", "owner3_share", "owner4_share", "owner5_share")
subset_ownership<-subset_ownership %>% 
  dplyr::mutate(own_mean = rowMeans(dplyr::select(., shares), na.rm=TRUE),
         owner_single=ifelse(owner1_share==100,1,0))

# Merge COMPUSTAT & ORBIS data
firms.de.ownership<-merge(firms.de.analysis, subset_ownership, by=c("conm","fyear"), all = TRUE)


# TABLE B.1 Models 1- 3
attach(firms.de.ownership)
lm_owner_1<-rdrobust(owner1_share, emp, c=2000)
lm_owner_1$coef["Bias-Corrected",]
lm_owner_1$ci["Bias-Corrected",]
lm_owner_1$pv["Bias-Corrected",]

lm_owner_mean<-rdrobust(own_mean, emp, c=2000)
lm_owner_mean$coef["Bias-Corrected",]
lm_owner_mean$ci["Bias-Corrected",]
lm_owner_mean$pv["Bias-Corrected",]

lm_owner_single<-rdrobust(owner_single, emp, c=2000)
lm_owner_single$coef["Bias-Corrected",]
lm_owner_single$ci["Bias-Corrected",]
lm_owner_single$pv["Bias-Corrected",]

# TABLE B.1 Models 4- 8
attach(firms.de.analysis)
lm_man<-rdrobust(man, emp, c=2000)
lm_man$coef["Bias-Corrected",]
lm_man$ci["Bias-Corrected",]
lm_man$pv["Bias-Corrected",]

lm_service<-rdrobust(service, emp, c=2000)
lm_service$coef["Bias-Corrected",]
lm_service$ci["Bias-Corrected",]
lm_service$pv["Bias-Corrected",]

lm_tech<-rdrobust(tech, emp, c=2000)
lm_tech$coef["Bias-Corrected",]
lm_tech$ci["Bias-Corrected",]
lm_tech$pv["Bias-Corrected",]

lm_trade<-rdrobust(trade, emp, c=2000)
lm_trade$coef["Bias-Corrected",]
lm_trade$ci["Bias-Corrected",]
lm_trade$pv["Bias-Corrected",]

lm_year<-rdrobust(as.numeric(as.character(fyear)), emp, c=2000)
lm_year$coef["Bias-Corrected",]
lm_year$ci["Bias-Corrected",]
lm_year$pv["Bias-Corrected",]


#------------------------------------------------------------------------------------------------------#
# 4) FIGURE B.7: MCCRARY DENSITY TEST PLOT
#------------------------------------------------------------------------------------------------------#

DCdensity(emp, 2000, bin = NULL, bw = NULL, verbose = FALSE, plot = F, ext.out = FALSE, htest = FALSE)

#------------------------------------------------------------------------------------------------------#
# 5) TABLE B.2: ROBUSTNESS EXCLUDING FIRMS WITH FOREIGN SUBSIDIARIES
#------------------------------------------------------------------------------------------------------#

fsubsidiary<-read_excel("~/Dropbox/02_PhD/07_Corporate_Savings/03Paper/05Submission_Draft/foreign_subsidiaries.xls") %>% 
  dplyr::select(-emp)
firms.subsidiary<-merge(fsubsidiary, firms.de.analysis, by =c("conm", "fyear"))
firms.de.analysis.nf<-subset(firms.subsidiary, foreign!=1)

attach(firms.de.analysis.nf)

# TABLE B.2: - MODEL (1) - Baseline
model1<-(rdrobust(cash, emp, c=2000, all=TRUE))
model1$coef["Bias-Corrected",]
model1$pv["Bias-Corrected",]
model1$ci["Bias-Corrected",]

# TABLE B.2:- MODEL (2) - Adding year/ indsutry-fixed effects
model2<-(rdrobust(cash, emp, c=2000, covs = cbind(as.factor(fyear), man, service, tech, trade)))
model2$coef["Bias-Corrected",]
model2$pv["Bias-Corrected",]
model2$ci["Bias-Corrected",]

# TABLE B.2: - MODEL (3) - Adding year/ indsutry-fixed effects
model3<-rdrobust(cash, emp, c=2000,
                 covs = cbind(as.factor(fyear), man, service, tech, trade),
                 cluster=gvkey)
model3$coef["Bias-Corrected",]
model3$pv["Bias-Corrected",]
model3$ci["Bias-Corrected",]

#------------------------------------------------------------------------------------------------------#
# 6) TABLE B.3: ROBUSTNESS FOR DIFFERENT MEASURES OF CORPORATE SAVINGS
#------------------------------------------------------------------------------------------------------#

attach(firms.de.analysis)

# TABLE B.2:- MODEL (1) - Logged savings 
model1<-rdrobust(ln_cash, emp, c=2000, covs = cbind(as.factor(fyear), man, tech, service, trade))
model1$coef["Bias-Corrected",]
model1$pv["Bias-Corrected",]
model1$ci["Bias-Corrected",]
model1$bws

# TABLE B.2: - MODEL (2) - Logged savings & clustered se
model2<-rdrobust(ln_cash, emp, c=2000,
                 covs = cbind(as.factor(fyear), man, tech, service, trade),
                 cluster=gvkey)
model2$coef["Bias-Corrected",]
model2$pv["Bias-Corrected",]
model2$ci["Bias-Corrected",]
model2$bws

# TABLE B.2: - MODEL (3) - Savings ratio
model3<-rdrobust(cash_ratio, emp, c=2000,
                 covs = cbind(as.factor(fyear), man, tech, service, trade))
model3$coef["Bias-Corrected",]
model3$pv["Bias-Corrected",]
model3$ci["Bias-Corrected",]
model3$bws

# TABLE B.2: - MODEL (4) - Savings ratio & clustered se
model4<-rdrobust(cash_ratio, emp, c=2000,
                 covs = cbind(as.factor(fyear), man, tech, service, trade),
                 cluster=gvkey)
model4$coef["Bias-Corrected",]
model4$pv["Bias-Corrected",]
model4$ci["Bias-Corrected",]
model4$bws


#------------------------------------------------------------------------------------------------------#
# 7) TABLE B.4: ROBUSTNESS FOR SAVINGS PER WORKERS
#------------------------------------------------------------------------------------------------------#

attach(firms.de.analysis)
# TABLE B.2: - MODEL (1) - Baseline
model1<-(rdrobust(cash_emp, emp, c=2000, all=TRUE))
model1$coef["Bias-Corrected",]
model1$pv["Bias-Corrected",]
model1$ci["Bias-Corrected",]
model1$bws

# TABLE B.2:- MODEL (2) - Adding year/ indsutry-fixed effects
model2<-(rdrobust(cash_emp, emp, c=2000, covs = cbind(as.factor(fyear), man, service, tech, trade)))
model2$coef["Bias-Corrected",]
model2$pv["Bias-Corrected",]
model2$ci["Bias-Corrected",]

# TABLE B.2: - MODEL (3) - Adding year/ indsutry-fixed effects
model3<-rdrobust(cash_emp, emp, c=2000,
                 covs = cbind(as.factor(fyear), man, service, tech, trade),
                 cluster=gvkey)
model3$coef["Bias-Corrected",]
model3$pv["Bias-Corrected",]
model3$ci["Bias-Corrected",]

#------------------------------------------------------------------------------------------------------#
# 8) FIGURE B.9: ALTERNATIVE RDD SPECIFICATIONS WITH CHANGING BANDWIDHTS
#------------------------------------------------------------------------------------------------------#

b <- 10:200
container <- data.frame()
container2<-data.frame()
bandwidths<-data.frame(rep(10:200, each=3))

attach(firms.de.analysis)

# Conduct RDD for different manual bandwidths
for(i in b){
  a <- rdrobust(cash, emp, c=2000, h=i)
  container <- rbind(container, a$coef)
  container2<-rbind(container2, a$se)
  sensitivity<-cbind(container, container2)
}

sensitivity<-cbind(sensitivity, bandwidths)
colnames(sensitivity)<-c("coef","se","bandwidth")
ind <- seq(1, nrow(sensitivity), by=3)
conventional_bw<-sensitivity[ind,] %>% mutate(upper=coef+1.6465*se,
                                              lower=coef-1.6465*se,
                                              upper95=coef+1.96*se,
                                              lower95=coef-1.96*se,
                                              optimal=ifelse(bandwidth==164,1,0))

# Display models as plot
ggplot(conventional_bw, aes(x = bandwidth, y = coef)) + geom_line(color="black", size=0.8) +
  geom_ribbon(aes(ymin=lower, ymax=upper), alpha=.5, fill="grey5") +
  geom_ribbon(aes(ymin=lower95, ymax=upper95), alpha=.4, fill="grey") +
  geom_hline(yintercept = 0, color= "black", linetype = "dashed") +
  theme(text=element_text(size=12, family="Frutiger-Light"), legend.position = "") +
  xlab("Bandwidth of local linear regression") + ylab("Estimated Effect") +
  geom_point(data=subset(conventional_bw, conventional_bw$optimal==1), shape = 21, fill = "black", color = "black", size = 1.5) +
  geom_errorbar(data=subset(conventional_bw, conventional_bw$optimal==1), aes(ymin = lower95, ymax = upper95), 
                height = .001, width = 1,  color = "grey3")

#------------------------------------------------------------------------------------------------------#
# 9) FIGURE B.9: FIRM-LEVEL DIFFERENCE-IN-DIFFERENCE TESTS
#------------------------------------------------------------------------------------------------------#

firms.panel<-firms.de %>% filter(fyear>1990 & fyear<2018) %>%
  mutate(per=ifelse(emp>=2,1,0), # generate treatment dummy
         cash=ch + ivst, # generate cash variable 
         emp=emp*1000, 
         sic_one=substr(sic,1,1)) %>% # create one-digit SIC code
  filter(sic_one!=6, !grepl(' SE', conm), # drop financial firms 
         !gvkey %in% c("236260","242158","216075")) # drop SEs


firms.panel<- firms.panel %>% group_by(conm) %>% dplyr::mutate(emp.beginning=first(emp)) %>% 
  filter(emp.beginning<2000 & emp.beginning>500) %>% # drop alwaystakers
  mutate(per.lead1=dplyr::lead(per, n=1L), # create leads, lags and changes
         per.lead2=dplyr::lead(per, n=2L),
         per.lead3=dplyr::lead(per, n=3L),
         per.lag1=dplyr::lag(per, n=1L),
         per.lag2=dplyr::lag(per, n=2L),
         per.lag3=dplyr::lag(per, n=3L),
         emp.lead1=dplyr::lead(emp, n=1L),
         emp.lead2=dplyr::lead(emp, n=2L),
         emp.lead3=dplyr::lead(emp, n=3L),
         emp.lag1=dplyr::lag(emp, n=1L),
         emp.lag2=dplyr::lag(emp, n=2L),
         emp.lag3=dplyr::lag(emp, n=3L),
         emp.lag4=dplyr::lag(emp, n=4L),
         emp.change1=emp-emp.lag1, 
         emp.change2=emp.change1-emp.lag2,
         emp.change3=emp.change2-emp.lag3,
         emp.change4=emp.change3-emp.lag4,
         emp.change.lead1=emp.lead1-emp,
         emp.change.lead2=emp.lead1-emp.lead2,
         emp.change.lead3=emp.lead2-emp.lead3)

terms<-c("per.lead3" , "per.lead2" , "per.lead1", "per", "per.lag1", 
         "per.lag2", "per.lag3")

adoption<-felm(cash ~  per.lead2 + per.lead1 +
                 per + per.lag1 + per.lag2 + per.lag3  + emp  + emp.change.lead2 + emp.change.lead1 +
                 emp.change1 + emp.change2 + emp.change3 + sale  | fyear  + conm  | 0 | conm, 
               data=firms.panel) %>% tidy %>% unnest %>%  filter(term %in% terms) %>% 
  mutate(lo=(std.error*qnorm(0.025))+estimate, hi=(std.error*qnorm(0.975))+estimate) 


ggplot(adoption, aes(x=factor(term), y=estimate)) + 
  geom_point(position=position_dodge(0.6)) +
  scale_x_discrete(limits=c("per.lead2","per.lead1","per", "per.lag1","per.lag2","per.lag3"),
                   labels=c("2 Years \n Prior", "1 Year \n Prior",
                            "Year of \n Adoption","1 Year \n After","2 Years \n After",
                            "3 Years \n After")) +
  geom_linerange(aes(ymin=lo, ymax=hi),position=position_dodge(0.6)) + 
  theme_bw() + geom_hline(yintercept=0, linetype="dashed") + xlab("") + ylab("Corporate Saving") +
  theme(text=element_text(size=12, family="Frutiger-Light"))




